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Abstract 

Two-dimensional (2-D) homogeneous magnetohydro dynamic (MHD) turbulence has many of 
the same qualitative features as three-dimensional (3-D) homogeneous MHD turbulence. These 
features include several ideal invariants, along with the phenomenon of broken ergodicity. Broken 
ergodicity appears when certain modes act like random variables with mean values that are large 
compared to their standard deviations, indicating a coherent structure or dynamo. Recently, the 
origin of broken ergodicity in 3-D MHD turbulence that is manifest in the lowest wavenumbers was 
explained. Here, a detailed description of the origins of broken ergodicity in 2-D MHD turbulence 
is presented. It will be seen that broken ergodicity in ideal 2-D MHD turbulence can be manifest in 
the lowest wavenumbers of a finite numerical model for certain initial conditions or in the highest 
wavenumbers for another set of initial conditions. The origins of broken ergodicity in ideal 2-D 
homogeneous MHD turbulence are found through an eigenanalysis of the covariance matrices of 
the modal probability density functions. It will also be shown that when the lowest wavenumber 
magnetic field becomes quasi-stationary, the higher wavenumber modes can propagate as Alfven 
waves on these almost static large-scale magnetic structures. 

PACS numbers: 05.20.Jj, 47.27.Gs, 91.25.Cw, 95.30.Qd 
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I. INTRODUCTION 


Turbulence in conducting fluids is a universal phenomenon occurring in stars, stellar 
winds, planetary interiors and magnetospheres, as well as in magnetic fusion energy devices 
and in the liquid metal cooling systems of fission reactors. Understanding such magnetohy- 
drodynamic (MHD) turbulence is of great importance because of the enhanced small-scale 
transport it affords when compared to molecular viscosity, resistivity and thermal conduction, 
in addition to the processes of large-scale self-organization it enables through its inherent 
nonlinearity. Although much progress has been made 1 , the study of turbulence is still a chal- 
lenging endeavor, dne to the vast number of independent variables that must be modeled or 
measured, which may be characterized as turbulent eddies and magnetic flux tubes, or as 
Fourier modes corresponding to length scales ranging from the size of the confining system 
to the dissipation scale at which nonlinear effects are no longer important. As is well known, 
the effective number of degrees-of-freedom in three-dimensional (3-D) flows is proportional 
to Reynolds number to the 9/4 power and in two-dimensions (2-D), it is proportional to 
Reynolds number squared 2 . Although this number is finite, it is large enough that numerical 
simulation of turbulent systems remains very challenging. 

Turbulence is thus a problem of statistical physics 3,4 and a firm basis from which to 
begin understanding MHD turbulence is equilibrium statistical mechanics 5,6 , which, by its 
nature, requires us to consider non-dissipative, conservative model systems. Although a gas 
of atoms or molecules may have only one essential conserved quantity (the total energy), 
in the case of ideal (that is, non-dissipative) homogeneous MHD turbulence there are three 
global invariants 5,6 (unless a constant external mean magnetic field is imposed, in which 
case it has two invariants 7 ). A gas, with only one ideal invariant, is expected to have this 
energy equipartioned over its degrees-of-freedom, while the presence of more than one ideal 
invariant in either 2-D or 3-D MHD turbulence leads to a statistical distribution of energy 
that may strongly favor a small fraction of available modes 5,6 . These favored modes are 
at the largest length scales in 3-D and may be at either the largest or the smallest length 
scales in models of ideal 2-D MHD turbulence, as will be seen. When these favored modes 
have sufficiently more energy than other modes, broken ergodicity (he., apparently non- 
ergodic behavior that does not disappear on reasonably long timescales 8 ) occurs and the 
magnetofluid ‘self-organizes’ so that these modes appear to behave like random variables 
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with large mean values and relatively small standard deviations 9 . In other words, what may 
be called a magnetic dynamo or coherent structure occurs naturally in ideal MHD turbulence 
due to the process of ‘broken ergodicity’. 

Broken ergodicity and its relation to magnetic dynamos and coherent structure was dis- 
cussed recently in great detail for 3-D homogeneous MHD turbulence, where eigenanalysis 
was used to discover the essence of the phenomena 10 . Here, we focus on ideal 2-D homo- 
geneous MHD turbulence and show that the 2-D case contains novel features not found in 
the 3-D case. We study a finite Fourier model of 2-D MHD turbulence and determine the 
eigenvariables and eigenvalues associated with the modal Hcrmitian 2x2 covariance matrices 
appearing in the probability density function. This allows us to analyze data from numer- 
ical experiments and uncover the underlying structure of 2-D ideal MHD turbulence, and 
also to observe the propagation of higher wavenumber Alfven waves on quasi-static lowest- 
wave number magnetic fields even when there is no mean magnetic held present. Also, 
although Fourier models have periodic boundaries, these can, serve as surrogates of physi- 
cal confinement 11-13 , allowing an investigation of large-scale behavior. The restriction to a 

2- D planar geometry both reflects the observed 2-D behavior of plasmas in relatively strong 
external magnetic fields 14,15 and allows access to much higher wavenumbers for long-time 
numerical simulations than are possible in 3-D simulations. In this paper, ideal 2-D MHD 
turbulence is of primary concern and although broken ergodicity in dissipative MHD turbu- 
lence can also be studied using 2-D models, which could be very informative even though 

3- D models are more realistic, this will be deferred for now. 

II. BASIC EQUATIONS 

Here, a turbulent magnetofluid will be assumed to how only in the x-y plane, with position 
vectors x = xx + y y, where x and y provide an orthonormal basis and where z is a unit 
vector normal to this plane, so that z • x x y = 1- The non-dimensional form of the 2-D 
incompressible MHD equations are well known 1 , and may be written as (with dt = d/d t, 
etc.), 

d t u: = Z (j, a) + Z(ip, oj) + B 0 ■ Vj + zA7 2 tn, (1) 

d t a = Z(ip, a) + B 0 • Vip + ?/V 2 a. (2) 
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The Jacobian Z(f,g ) contains the nonlinear interactions and has the form 

Z(f, g ) = d x f d y g - d y f d x g. (3) 

The vorticity u;(x, t) is related to the stream function ip(x, t) by uj = —W 2 ip, and the electric 
current j(x, t) is related to the magnetic potential a(x, t) by j = — V 2 a. The magnetofluid 
velocity is u = V x zip, the magnetic induction is b = V x z a, and these clearly satisfy 
V ■ u = V • b = 0, where V = x<9 x + yd y . Density does not appear in (1) because it equals 
unity, v in (1) is the kinematic viscosity, while rj in (2) is the magnetic diffusivity. The 
variables ip(x,t), a(x, t) and their derivatives are assumed to be periodic within a square 
of edge length 27r; this periodic square can be used to represent a small region within a 
much larger area of homogeneous flow or as a surrogate for a bounded physical area (or as 
a manifold without boundary). In eqs. 1 and 2, B 0 is the constant, mean magnetic field 
(which can be zero). 

Here, we assume periodic boundary conditions on a square with side length 27 r. Let us 
define the spatial average of a quantity /(x)g(x) as fg, where 

r 2n r 2n 

fg=^ 2 j o dy Jv dxf(x.)g(x). (4) 

Spatial averages of the energy E, cross helieity He and mean squared magnetic potential A, 
as well as the enstrophy and mean square current J, are defined by 

E = ±(ipu + aj) = !(w 2 + b 2 ), (5) 

H c = ioa) = iu • b = i jip, (6) 

A - A n = p, J = if. (7) 

We define these averages first, rather than the total quantities, because these averages are 
independent of numerical grid size. Initial conditions are generally chosen so that E — 1, 
as is appropriate for the dimensionless equations (1) and (2). On an N x N numerical 
grid, the total energy is E = N 2 E, and similarly He = N 2 He, A = N 2 A, etc. Spatial 
grid-point averages are useful in that they allow a more straightforward comparison between 
simulations with different values of N. 
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Using equations (1) and (2), it is straightforward to find 


§ = - 2 (^ + V J\ 


( 8 ) 


dH, 


c 


dt 


[v + rj)ju, 


(9) 


dA 

dt 


= B 0 x z • du — r/b 2 . 


( 10 ) 


ff we have v — rj — 0, then E and He are constants of the motion, and if B 0 = 0, then A 
is also constant. In de-aliased Fourier models of 2-D MHD turbulence, it can be shown that 
these three quadratic forms are the only constants of the motion 16 . 

Briefly, consider the case where B 0 = 0 and r) = 0, and the equation (2) takes the form 


da da „ 

_ = — + u • Va = 0. 


( 11 ) 


This indicates that the value of a for every moving fluid point in an ideal 2-D magnetofluid 
is a constant of the motion. However, to capture this property, we must adopt a Lagrangian 
approach and use a numerical method that tracks individual points. Instead, the standard 
procedure in homogeneous MHD turbulence simulations is to follow an Eulerian approach 
that uses a Fourier spectral transform method 17 , a method based on the time-dependent 
values of a (and co) on a fixed, regularly spaced, finite grid of points. This Eulerian approach 
does not allow the tracking of individual fluid elements, so, in general, the Lagrangian 
invariants a(x(t)) associated with each fluid element play no role in the Fourier method. 
Although a Lagrangian approach that tracks individual vortex lines exists for ideal 2-D fluid 
mechanics 18 , where the governing equation is (11) with a replaced by co (he., dco/dt = 0), 
and for the analogous 2-D motion of parallel electric line charges 19 , there appears to be no 
similar approach possible in ideal 2-D MHD, as this would require du/dt = 0 and dj /dt = 0, 
and neither of these are the equations of 2-D MHD. 


III. FOURIER MODELS 

The system we employ to model 2-D MHD turbulence uses the finite Fourier series cor- 
responding to co and a. The physical fields ca(x) and a(x) in x-space are connected to their 
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corresponding Fourier coefficients u)(k) and a(k) in k-space by transformations of the form 
(where i = \/— T) 

= (12) 
ke/c xe(P 

The sums are over discrete sets /C and V. The set /C contains those wave vectors (he,, k-space 
grid points) k that have integer components k x = k\ and k y = k 2 , whose values he between 
— N/2 + 1 and N/ 2; /C thus has N 2 elements. The set V consists of position vectors (he., 
x-space grid points) x that have components x 3 = 2ix ntij/N (j — 1,2; 0 < rrij < N ), with 
x = Xi , y = x 2 , rn x = rn \ , and rn y = m 2 , where the rn 3 are integers; N is the number of grid 
points in each dimension and the set V also has N 2 elements. 

Since x-space variables, such as a(x) are real, their Fourier transforms satisfy relations of 
the form a(— k) = a*(k), where ‘*’ indicates complex conjugation. The a(k) have real and 
imaginary parts a^(k) and a/(k), so that a(k) = a^(k) + iaj(k), which implies 

flfl(-k) = a R { k), a/(-k) = -a 7 (k). (13) 

Thus, only about half of the a(k) for k G /C are independent, i.e., those corresponding 
to k, but not to — k. Now, each a(k), k 6 1C, generally has one real and one imaginary 
component, i.e., except for coefficients with k = (k x ,k y ) = (0,0), (0, iV/2), (N/ 2,0) and 
(N/2, N/2), which only have real components. Nevertheless, a careful count shows that 
the set of coefficients a(k), k e /C, also have a N 2 independent parts, before any further 
restrictions are imposed. 

However, these excepted coefficients are not used as dynamical variables in Fourier method 
solutions of equations (1) and (2), because only Fourier coefficients with k = |k| such that 
1 < k < K < N/2 are allowed to be nonzero. Coefficients with wavevectors k outside 
the ‘isotropic truncation radius’ K , or at k — 0, will always be set to zero. (Here, we use 
K = y/2N/Z, a value set by the de-aliasing procedure 20 we employ.) The largest dynamical 
scale in the periodic box is represented by the k — 1 modes and although /C contains a 
total of N 2 wavevectors, only A f = nK 2 vectors k satisfy 1 < k < K and participate in the 
dynamical evolution of the model system. For convenience, let us define the set 1C C /C by 

K! = {k | 1 < k < K; if k G 1C, then - k £ K'} . (14) 

The set 1C has A f = J\f / 2 elements. If we have a set 1C, an equivalent set can be created 
by removing k from 1C and replacing it by — k; there are obviously 2 A/7 possible equivalent 
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choices for K! . The exact choice depends on the particular computer code we use; in what 
follows, for k = (k x , 0), we allow 1 < k x < K\ and for k = ( k x ,k y ), we allow 1 < k y < K 
and —A' <k x < K, so long as 1 < ( k 2 + ky) < K. 

Using (12), it can easily be shown that for any pair of variables / and g, 


E/OmMm) = E/*(M)2(M)- 


x£(P k &K 

In addition, the global average (4) can now be expressed as 


(15) 


fa = /‘(k,«)g(k,«). (16) 

k etc 

In what follows, we often omit t from the argument of the coefficients. 

Performing a Fourier transform on the equations (1) and (2), yields a finite, autonomous 
dynamical system: 


dto( k) 
dt 

da( k) 
dt 


{Z(j.“)}k + |Z!t', J)} k - >B., ■ k/iki - vk 2 d i(k), 
{Z(ip,a)} k + 1B 0 - ki/>(k) - t)k 2 a( k). 


( 17 ) 

(18) 


Here, we use (12) to transform the relations to = — V 2, 0 and j = — V 2 a into k-space: 


o)(k) = k 2, ip( k), j(k) = /c 2 a(k). 


(19) 


Terms such as {Z(j,a)} k denote the Fourier transform of Z(j,a), etc. 

The dynamical system defined by (IT) and (18) contains the quadratic nonlinear terms 
{Z(j, a)} k , {Z('0, ca)} k and {Z^, a)} k . This dynamical system is high-dimensional for 
J\f » 1 and numerical solutions for small values of v and g are highly chaotic, i.e., turbu- 
lent. An efficient algorithm for numerical integration of (17) and (18) is the Fourier spectral 
transform method 17 coupled with a third-order time-integration scheme 21 . The essence of 
this algorithm is to integrate the equations (IT) and (18) in k-space and to determine the non- 
linear terms such as {Z(j, a)} k , when needed, by Fourier transforming the separate quadratic 
factors back to x-space, forming products, and then transforming forward to k-space. For ex- 
ample, the first factor of Z(f, g ) appearing in (3) is found from /(k) using the first transform 
in (12): 

s,/ = ijiy(k)t' t " (20) 

ke£ 
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Once all the factors are transformed and Z(f,g) is determined at each x-space grid point, 
we transform back to k-space: 

{Z(/,0)} k = (21) 

x6(P 

Note that when the transform (21) is performed, all coefficients {Z(/, <?)} k outside the range 
1 < k < K are immediately set to zero (where K = \/2N/3). The transforms (20) and (21) 
are done using fast Fourier transform (FFT) routines for speed, and the whole procedure 
is actually done twice, once to the x-space grid and once to a shifted x-space grid, with 
the results being averaged to remove any remaining aliasing error 20 . This Fourier spectral 
transform method thereby confines the model dynamical system to the finite set of Fourier 
coefficients u)(k) and a(k) for which 1 < k < K. The algorithm would produce exact 
solutions except for time-integration and round-off errors that introduce random, quasi- 
thermal fluctuations into the model dynamical system. 

These fluctuations, in fact, serve as a ‘heat bath’ so that the statistical mechanics of 
our computer model can be described in terms of canonical ensembles. In the next section 
we discuss the statistical theory of ideal 2-D MHD turbulence. Following this, we present 
numerical examples and then elucidate the origin of broken ergodicity in the dynamical 
system (17) and (18) that models 2-D MHD turbulence. Finally, we discuss our results and 
conclude this paper. 

IV. STATISTICAL MECHANICS 

Statistical studies of MHD turbulence represented by Fourier modes were initiated by T. 
D. Lee in an early paper 22 that demonstrated the existence of canonical ensembles based on 
ideal invariants, as well as a Liouville theorem. The possible constants of the motion for 
2-D MHD are E, He and A, as defined in (5), (6) and (7), respectively. The term ‘ideal 
invariants’ will be used for the quantities E = N 2 E, He = N 2 Hc and A = N 2 A; using (16), 
these take the forms 

E = k), -F(k) = k~ 2 \uj(k)\ 2 + /c 2 |a(k)| 2 , (22) 

k eic 

H c = J2 H c( k ), H c Qs.) = u R (k,t)a R (k) + u)j(k, i)aj(k), 

k GK' 
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(23) 



( 24 ) 


A = 2>(k), 2l(k) = |S(k)| 2 . 

k e/C' 

Here, the summations are over the set of k dehned in (14), and E, He and A are the total, 
rather than the average energy, cross helicity and squared magnetic potential; similarly, total 
enstrophy is = N 2 Vt and total squared current is J = N 2 J. The quantities E, Hq and A 
are used in the ‘absolute equilibrium ensemble theory’ of 2-D MHD 6 . 

This statistical mechanics of ideal 2-D MHD turbulence is based on a Gaussian canonical 
probability density function (PDF), D, which can be represented as 

D = — exp(— aE — /3Hc — 7 A) (25) 

Zj 

Z = J exp(— ck E — j3Hc — , yA)dr, (26) 


dr = n duji{k)dui(k.)daR(k.)ddj(k). (27) 

ke/C' 

In the equations above, Z is the total partition function and dT is the phase space volume 
element dehned by the independent variables contained in the Fourier model (he., k e JC); 
each independent variable is integrated from — 00 to + 00 . The set KJ contains A f = \ix K 2 
wave vectors k that satisfy 1 < k < K, and the associated independent variables are a)^(k), 
u)j(k), o,r( k) and d/(k), for a total of Nr = 4J\f = 2 ttN 2 coordinates in the phase space T. 
The PDF also applies to the case where B 0 ^ 0, if we set 7 = 0 in (25) and (26). 

The PDF (25) can be used to produce ensemble predictions (<F) of phase functions <h, 
which may be compared with time- averages <h: 


($) = f $Ddr, $ = 1 f T $dt. 

Jr T Jo 


(28) 


The phase functions in (28) are moments or sums of moments of these Nr phase variables. 
Ergodicity is informally dehned as occurring when (<3>) = <h for all <h, and non-ergodicity as 
occurring when ($) ^ $ for some $. 

The quantities D and Z in (25) and (26) can be expressed in terms of the modal quantities 
c5(k) and a(k) as follows: 

ex p[— <3(k)] 


D = n r>(k). n(k) = 

ke/C' 


Z(k) 


( 29 ) 
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(30) 


Z 


II Z( k), Z(k) = / exp[-Q(k)]dr(k), 

k &K.' J 


t/r (k) = duR(k)dui(k)daR(k)dai(k). (31) 

Again, the product in (29) is only over independent k and all integrations in (30) are from 
— oo to +oo. The series (22), (23) and (24) are used in (25), and thence (29), to produce 
Q(k), a necessarily positive definite, real quadratic form 

Q(k) = aE{ k) + f3H c ( k) + 7 A( k). (32) 


where the expressions for E( k), Hq{ k) and A(k) dehned in (22), (23) and (24) can be placed 
into Q(k), and the modal PDF (30) can then be expressed in terms of u)(k) and a(k) as 

D{ k) = ^ _1 (k)exp[-y t (k)M fc jjr(k)]. (33) 


Here, ?/'(k) is the Hermitian adjoint of the column vector y( k), where 

y ( k ) = 


£l( k ) 


uj(k)/k 

1/2 (k) 


k a(k) 


The real, symmetric (he., Hermitian) 2x2 covariance matrix M k is 

0/2 


Mi, = 


a 


P/2 a + ^/k 2 

Now, our next step is to hnd the eigenvalues and eigenvectors of M k . 


(34) 


(35) 


A. Eigenanalysis 


Diagonalization of the modal matrices M k in (35) is straightforward and proceeds via a 
similarity transformation using a unitary matrix T k , where 


T, = 


G k P 
-P G k 


y/2C k G k 

Application of T k yields the diagonal matrix A k : 
Afc = T k M k T k = 


. = + r G t = c k ~p r 


1 

V 

O 

i 


a + h (f + Gk) 

0 

i 

O 

V 

1 


0 

a + h (*r - C k ) _ 


(36) 


(37) 
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Here, T k is the Hermitian adjoint of T k , and Aj. 1,2 '* are the eigenvalues of M k . Although T k is 
really orthogonal and T\ is simply its transpose, the vectors ?/(k) and eigenvectors h(k) are 
complex, we can class T k as unitary. The eigenvectors h(k) are 


h(k) 




{^(k) 

1 

G k Cj{k)/k + j3ha{ k) 

tA 2 )(k) 

p2C k Gk 

—(3Cj(k)/k + G k ka( k) 


(38) 


The complex quantities tA 1,2 )(k) will be called eigenvariables. Note that if B 0 p 0, then 
7 A = 0 in (25) and the tA 1,2 )(k) become the Elsasser variables 23 , ^ ± (k), in 2-D MHD: 


lim {;d> 2 )(k) 
7.4— »o 



± Kfk, 
k 


1 

71 


w(k) 

k 


=F fca(k) 


±^^ ±(k) . 


P> 0, 



p< 0. 


(39) 


Thus, the eigenvariables fA 1,2 )(k) encompass the ^ ± (k) as a special case when B 0 p 0. 

In terms of the real and imaginary components of the eigenvariables, Vr’ 2 ^ (k) and hj 1,2) (k), 
the modal PDF (30) takes the form 

D(k) = D (1) (k)D (2) (k), D (1 > 2) (k) = D2’ 2) (k)D?’ 2) (k), (40) 


Z(k) = Z (1) (k)Z (2) (k), Z (1 > 2) (k) = zg’ 2) (k)zj 1,2) (k), (41) 


4 1,2, (k) = 




41,2) 


(k) 


S = R, I, 


(42) 


4 1,2) ( k ) = 


/ exp(-A7 2) [4 1 ’ 2) (k)] 2 )df4 1 ’ 2) (k) = /-^ y 


(43) 


We can use these results to find expectation values of the modal eigenvariables, as well as 
expectation values for modal and total E, Hq : A, and other quantities, as needed. These 
expectation values will be rational functions of the ‘inverse temperatures’ a, (3 and 7, which 
are, at this point, undetermined. This, in turn, will lead to a relationship between a, (3, 
7 and E, He , A, a relationship that will be parameterized by a single real variable. A 
minimization procedure will then be introduced, one that gives a unique value to the single 
real variable and thereby determines a, [3 and 7, as well as the various modal expectation 
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values. Following this, we can compare these canonical ensemble predictions with time- 
averages from numerical simulations. There will be a disparity between the two, a broken 
ergodicity that can be completely understood using eigenvalues and eigenvariables. 


B. Expectation Values 

Now that the total PDF (25) has been expressed as the product of the PDFs (42) of the 
real and imaginary components of the eigenvariables (which provide coordinates for the phase 
space T), we can easily determine expectation values for moments of these eigenvariables. 
Using (42) and (43), it is easily shown that 

(®s' 2| (k)) = 0, («H’ 2, ( k )4 1 ' 2| ( k )) =0, (44) 


[U 2) ( k )] 2 } = 


2 A‘ WI 


, S = R,I. 


Using (22), (23), (24), (38), (44) and (45), we find 


<£(k)> = (|« ( 1 ) (k)| 2 ) + (|«< 2 >(k)| 2 ) = -i y + -E 

An An 


2a + 7 /k 2 

w~ 


(Hc( k)> = 


JL 

2C, 


z ' (1) ( k )| 2 ) — (p <2, ( k ) 


2 C k 


1 




(1) y(2) 


A 

‘2«r 


(MV) = 


c k - 7 A 2 
2 k 2 Ck 


h (1) (k)A + 


Ck + 7 A 2 
2 k 2 C k 


h (2) (k)p 


C k ~ 7 A 2 C k + 7//C 2 a 


2k 2 C k X ( k ) 2k 2 Ck\ ( ? ) k26 l ’ 


k^k 


5 ^ = a 2 — p /4 + « 7 //r 


In summary, we have 


(E(k)) = - \: /k \ (h c ( k» 


% 


OL 

2 J| ’ * * ' 


(45) 


(46) 


(47) 


(48) 

(49) 

( 50 ) 
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These are individual mode contributions to the expectation values of the ideal invariants. 
Note that, since He is a pseudoscalar under the parity transformation x — > — x, so is (3. 

We can also find the modal expectation value of quantities that are not ideal invariants, 
such as the magnetic energy; here 5% is given in (49): 

<£ M (k)> = k 2 (A(k))= (51) 

Note that > 0 — > = §1 > 0, and since (A(k)) > 0, we have a > 0. However, 

(50) tells us that /3 (H c { k)) < 0, so /3 takes a sign opposite to that of (H c ( k)). Also, since 
(m) > 0 for all k , then 7 > —2a, i.e., 7 can be positive, zero or negative. 

Now that we have various modal expectation values given in terms of a, [3 and 7, we 
can sum these and invert the results to find the heretofore unknown a, f3 and 7 in terms 
of the three integral invariants and one unknown parameter. We then use a minimization 
procedure to determine this unknown parameter and thereby determine a, [3 and 7, allowing 
us to assign values to ensemble predictions and to compare these to time averages found 
through numerical simulation. This, in turn, leads us to the mechanism of broken ergodicity. 

C. Inverse Temperatures 

Using (50) and (51), it is easy to verify that 


a(E(k)) + (3(H c (k))+'y(A(k))=2, 

(52) 

4 a 2 


a (E{k)) + — (H c (k))-’y(A(k))=0, 

(53) 

2a (HcQs.)) + p (E M (k)) = 0. 

(54) 


In what follows, we can set the expectation of values of the ideal invariants equal to their 
initial values, £, He and A, since these are expected to be constant, and set the expectation 
value of the magnetic energy equal to the variable <p, which is initially unknown, i.e., 

(E)=£, ( H C ) = H C , (A) — A, (E M ) = <p. (55) 

Now, we sum (52), (53) and (54) over the A f wavevectors k € 1C, using (55), to get 

a 2 

a£ + (3Hc + 7-4 = 2J\f, a£ + A— He — 7-4 = 0, 2 aHc + (3<p = 0. (56) 

r' 
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The equations (56) can be readily solved to yield the inverse temperatures in terms of the 
quantities £, He, A and p\ in terms of the grid-averaged values, e.g., p = p/N 2 , we have 


a = 


nip 


p(£ -p) - He 


2 > 


Uc £ — 2p 

p = — 2 ^^ck, 7 = — — a. 

P A 


(57) 


Here, n = J\f'/N 2 , where J\T' = nK 2 /2 and K = \/2N/3, so that n = 7t/9 = 0.3491. Thus, 
the undetermined multipliers a, [3 and 7 found in the PDF (25) are, in fact, all functions of 
one initially unknown variable p. Now, \Hc\ < £ and 0 < A < (p < 1 , and if £ — 1 , then a, 
f3 and 7 are all roughly of order unity, unless the denominators in the expressions are close 
to zero, which may occur. 


D. The Entropy Functional 


The parameters ck, (3 and 7 can be determined from the constant values S, T-Lc and A, 
as described by the relations (57), by minimizing an entropy functional 24,25 with respect to 
p. Here, we minimize the specific entropy functional cr(p), using the various formulae given 
above, 


& (<P) = ~w ( lnD ) = 2n ( 1 + ln7r ) - 7^2 ln5 fc ( 58 ) 

ke/C' 


5 


2 

k 


a 



£-2p\ 

Ak 2 ) ' 


(59) 


The form of <5 2 comes from using (49), along with the a(p), f3(p), p(p) and n, as defined 
in (57). Please note that, as discussed by Khinchin 24 , the entropy of a closed system, such 
as the finite Fourier model we are studying, is not a dynamic quantity that depends on the 
time-dependent value of p on a phase trajectory, but rather a global property of the ensemble 
phase space with a value of s = <j(p 0 ) determined by finding the p = p Q that minimizes 
cr(p). Using the fact that a > 0 and > 0, it can be shown that —00 < s < 1.7482 + 2 ln£. 

In an attempt to determine the minimum of cr(p), we can try to find the zeroes of its first 
derivative with respect to p. This first derivative is 

Idd = </(£) = F(0)G(0), (60) 

ap 
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(61) 


F{0) = 


if 3 - h c 2 ( 3 ^ - £J 

5 2 (p (S - ipj - He 2 


G (<p) = Jf2 E 


ip — Ak 2 


k e/C' £ — 2(p + Ak 2 1 - He /<p 2 


(62) 


Now, (u ■ b ) 2 < ( u 2 ) (6 2 ) implies He < <p{£ — 9?), so that the denominator of F(<p) must 

since /lu ± b| 2 \ > 0 


be positive. Then, (p_ < (p < (p +1 where <p± = \ ^ ± y £ 2 — AHc 

implies 2|^fc| < £, the (pA are real. The numerator of F (Cp) has a mimimum at (p = [Hc\ = 
yjip+ip- , at which value F {\pHc^ = He > 0, so that F(tp) > 0 for acceptable values of (p. 
Thus, the zeroes of cr'((p) are the zeroes of G{(p). 

Although there is no room here for an extensive discussion, there is only one value (p = <p 0 
such that G((p 0 ) = 0. Finding (p 0 requires a numerical procedure, except when B 0 ^ 0 (see 
below). Although G'(p) = dG/dtt > 0 can alternate in sign as (p is varied, we always have 
G\(p) > 0 at (p = ip 0 , where (after some algebraic manipulation), 


G\(po) = g'ipo) = Y, “ 


a k (1 + 3 h 2 /p 2 0 ) + 2 a 2 k h 2 /p\ 


ke/C' I 1 - 2 Po + H (1 - h 2 /p 2 0 )Y 

Po = Po/£, a k = Ak 2 j £ , h = \H c \/£. 


(63) 


The sum above contains terms that have a quadratic function of a k in the numerator. It 
can be shown that the quadratic is positive for all values of a k except possibly for a limited, 
inconsequential range of a k , while the corresponding denominators are always positive. The 
net result is that we always have g\p 0 ) > 0 and that (p = (p D minimizes cr(0), yielding the 
specific entropy s = cr(<p 0 ). Once (p Q is found, the inverse temperatures (57) are determined 
and ensemble predictions such as (50) can be assigned definite values. 


E. Nonzero Mean Magnetic Field 


The statistical mechanics of the case where B 0 ^ 0 has a much simpler formulation than 
the B 0 = 0 case given above. In the B 0 7^ 0 case, 7 = 0, or equivalently, (p Q = £/ 2, which 
leads to eigenvalues independent of k : 


A 1 - 2 ) 

A k 



2 n 


ST? H c 


(64) 
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In this case, the specific entropy is determined only by He (where, e = 2.718 . . . and n = tt/9): 

s = nln (In) - m c) • (65) 

Above, He < S/ 2, since we cannot have He — S/2 because then either z + (k) = 0 or 
z _ (k) = 0 for all k G /C', all nonlinear interactions are absent, and there is no turbulence. 
It is clear from the form of s in (65) that the maximum that s can take occurs at He = 0, 
where (again) s max = 1.7482 + 21n£, while the minimum s m i n is unbounded from below as 
He S/2. These limits are general for ideal 2-D MHD turbulence, as they also apply to 
the case where B 0 = 0, as already mentioned. 

There are n eigenvariables tA^(k) and n eigenvariables 7^ 2 ^(k) and together they have an 
average energy of E/(2J\f') = S/(2n). Individually, from (45), the eigenvariables have the 
following expectation values for their energy: 

(l« (1 ’ 2) ( k )l 2 ) = T<i3r < 66 > 

Afc 

In the B 0 ^ 0 and He = 0 case, s = s max and (65) and tells us that = A = S/(2n ), so 
that all eigenvariables have identical expected energies of S/{2n). However, if He — > S/2, 
then 4P — > oo, as seen in (64), and we have 

(^ (1) (k)) -G 0, (E (2) (k)) -> (67) 

Thus, when B 0 ^ 0 and He S/2, the Elsasser variables separate into one half with 
energies that are very large when compared to the other half. In the formalism given here, 
fA 2 )(k) is the half that always contains larger energies. Using (39) and (50), we see that for 
j3 > 0, the fA 2 )(k) are (except for a minus sign) the Elsasser variables z~( k) corresponding 
to He < 0, while for f3 < 0, the h^(k) are the £ + (k) corresponding to He > 0. 

V. NUMERICAL RESULTS 

As already mentioned, we numerically simulated the dynamical system defined by (17) 
and (18) by using a Fourier spectral transform method 17 coupled with a third-order time- 
integration scheme 21 . Here, we present twelve ideal runs (u — r) — 0) whose parameters are 
given in Table I. For all runs, B 0 = 0, N = 128, the time step size was At = 10~ 3 , and each 

run lasted 10 6 time steps, that is, from simulation time t — 0 to t — 1000. 
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TABLE I: Numerical Simulation Parameters. 


Run 

£ 

He 

A 

a 

P 

7 

<P 0 

s 

A 

1.00007 

0.35501 

0.42000 

2.2032 

-2.4189 

-1.5390 

0.64672 

1.0730 

B 

1.00000 

0.47309 

0.18632 

6.9125 

-12.717 

-1.0633 

0.51433 

0.80260 

C 

1.00006 

0.15709 

0.18632 

0.93966 

-0.50343 

-0.87158 

0.58644 

1.5687 

D 

1.00003 

-0.034993 

0.012160 

0.71012 

0.098225 

-0.69560 

0.50597 

1.7395 

E 

1.00000 

0.49931 

0.18632 

253.49 

-505.75 

-1.1105 

0.50041 

-0.46998 

F 

1.00006 

-0.0069888 

0.18632 

0.85806 

0.020222 

-0.85729 

0.59310 

1.6061 

G 

1.00005 

0.25147 

0.18632 

1.1090 

-0.96946 

-0.89648 

0.57533 

1.5010 

H 

1.00000 

0.47138 

0.00031952 

6.1330 

-11.953 

624.80 

0.48372 

0.92079 

I 

1.00000 

0.39998 

0.00015964 

1.8322 

-4.2739 

3605.3 

0.34293 

1.0983 

J 

1.00000 

0.38998 

0.00005499 

17.250 

-69.565 

192350. 

0.19340 

-0.87838 

K 

1.00000 

0.45960 

0.000091651 

20.207 

-58.490 

80445. 

0.31757 

-0.62948 

L 

1.00000 

0.38998 

0.00021996 

1.5870 

-3.1358 

1519.0 

0.39474 

1.2773 


In Table I, the values (to five significant figures) given for S, HLc and A are the time- 
averaged values over each complete run; at the beginning of each run, 8 = 1, and the time- 
averaged values differ by less than 1 part in 10 4 , indicating the accuracy of the numerical 
method (HLc an d A were similarly well conserved). These values of 8, HLc and A were used 
in (58) to find a, (3 and 7 by determining the value (p Q that minimized cr(0) to find the 
entropy s = cr((p 0 )] all these values, for each run, appear in Table I. 

Although 10 6 time steps might seem to be sufficient for the model dynamical system to 
reach equilibrium, Fig. 1 indicates that this is not always the case. In this figure, the values 
of mean square current J for Runs B, C, E, F and G (which all have A = 0.18632) are 
presented. While Runs C, F and G appear to have reached stationarity relatively quickly, 
and Run B soon thereafter, the value of J for Run E is still increasing. In looking at Table 
I, we see that Run E has the highest value of HLc (= 0.49931), which indicates that u«b, 
effectively bottlenecking the nonlinear transfer of energy from low to high wave numbers. In 
contrast, Run F has HLc — —0.007 and reaches its terminal value the quickest of all. 

In regard to the Runs B, C, E, F and G (all with A = 0.18632), Table I gives their 
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FIG. 1: Mean square current J for Runs B, C, E, F and G. 


respective values of He] correlating these values with the behavior shown in Fig. 1 indicates 
that the time needed to reach the terminal value of J is directly correlated to the magnitude 
of Tic- However, for Runs H, I, J, K and L, which all have relatively high values of [Hc\ but 
very low values of A, graphs analogous to Fig. 1 show these runs reached equilibrium quicker 
than any runs in Fig. 1. The remaining runs of Table I, Runs A and D, were similarly seen 
to have achieved equilibrium at about the same rate as Run G, shown in Fig. 1. 

A major qualitative difference of various runs in Table I is due to the value of Cp Q : All 
runs with (p Q > ± (7 < 0), such as Run A, peak at k — 1, while runs with (p 0 < \ (7 > 0), 
such as Run J, peak at k = K, the largest wavenumber in the simulation. The total 
ideal energy spectra for the eigenmodes can be calculated from the eigenvalues (37) and is 
E^ 2 \k) = nk/X^’ 2 ^-, example spectra for Runs A and J are shown in Fig. 2. Run A (c p Q > §, 
7 < 0) peaks at lowest k and Run J (<J5 0 < ±, 7 > 0) peaks at largest k ; if (p Q = \ exactly 
(7 = 0 ), as is the case when B 0 7 ^ 0 , the predicted spectra are precisely flat; if ~ §, as in 
Run E, then the expectation may be that the E^ l ' 2 \k) will be relatively flat, but it is not. 
Instead, the spectra for Run E are very similar to the spectra for Run A, as shown in Fig. 
2. Let us now take a more detailed look at the energy spectra for some of these runs. 

Two-dimensional spectra can be generated by determining the time-averages and vari- 
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FIG. 2: Ideal energy spectra for Runs A and J. 


ances, over t — 0 to 1000, for all of the Fourier coefficients of any given run. Average values 


and variances with respect to time will be denoted by u^(k) and |u^(k)|y (s = 1, 2), 
respectively; these are used to define the two-dimensional energy spectra 


EMk) = Ekftk)', £V(k) = £ 


/*>(k) 


S=1 


S=1 


2 

v’ 


E T (k) = E A (k) + E v (k). (68) 


These spectra are shown in Fig. 3 for Run E and are representative of the Cp Q > \ runs in 
Table I. Runs with (p Q > \ (such as Run J in Fig. 2) peak at k — K but have a similar 
relation between mean, variance and total two-dimensional spectra. 

The critical feature appearing in Fig. 3 is that E A { k) is not close to zero everywhere, but 
instead is relatively large at low k. In fact, at low k it gives the principal contribution to 
Et{ k). The statistical theory does not predict this phenomenon and something obviously is 
occurring that ‘breaks the ergodicity’ of the dynamical system under consideration. 

Returning to Fig. 2, what we see, for Runs A and J, is a comparison of the computed 
values (at t = 1000) of E^ l ' 2 \k) = nk\v^ 1,2 \k)\ 2 vg (here, ‘avg’ signifies values with the same 
k are averaged at a specific time t) with expected values (the dashed curves) determined 
by (66). The fit seems fairly good, particularly as k increases, and is representative of the 
level of closeness between predicted and computed values of F , ^ 1,2 '(/c) for the runs in Table I. 
However, at low k, there appears to be less of a match than as k — y K, and Fig. 3 indicates 
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FIG. 3: Run E, 2-D spectra. Energy in a. time-averages of the Fourier coefficients; b. variances 
(with respect to a.); c. total (sum of a. plus b.). 


that these unexpected values are not, in fact, momentary temporal fluctuation. We consider 
this further by now discussing the behavior if individual eigenmodes. 

In addition to recording instantaneous spectra, such as in Fig. 2, or examining time- 
statistics of spectra, as in Fig. 3, we also followed the time evolution of certain u^l(k) 
(s = 1, 2) for selected values of k. These can be used to produce time histories of E^ 1,2 \k), ; 
as well as phase trajectories, for these selected k. In Fig. 4, we present E^ 1,2 \k) vs time 
for k — 1, 2 and 4, from Run A; in this figure, we also have expected values for each of 
the E ( ' 1,2 \k). What we see is that, although the computed and expected E^ l ’ 2 \k) for k = 1 
fall on top of one another, those for k = 2 and 4 do not. Instead, the simulation values 
are anomalously high, commensurate with what we see in Fig. 2 at low k. Thus, although 
computed spectra appear much as predicted over moderate to large k , this is another clear 
indication that some form of non-ergodicity is present. 

In Figs. 5 and 6, both pertaining to Run A, we see more detail concerning the time 
histories shown in Fig. 4. These figures are ‘phase trajectories’, i.e., projections of the 
path the dynamical system is following in phase space onto a two-dimensional plane whose 
abscissa and ordinate are the real part and imaginary parts of a particular eigenvariable 
u( 1,2 )(k). Although Fig. 4 indicates that computed and expected values E^ 2 \k) for k — 1 
match, when we look the behavior of the individual coefficients in Fig. 5, we see that the 
mean magnitude of t^ 2 ^(k), k = (0,1) is significantly larger than expected (86.1 vs 66.9), 
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FIG. 4: Time histories of E^ s \k), s = 1, 2, for k = 1, 2 and 4, from Run A. 


while the mean magnitude of u^(k), k = (1,0) is significantly smaller than expected (37.0 
vs 66.9); the coefficients iA0(k) have close to the expected magnitudes (~ 0.59). Also, note 
that ?/ 2 )(k), k = (1,0) behaves like an Alfven wave, propagating on the almost stationary 
magnetic field associated with u®(k), k = (0, 1). 

In Fig. 6, we have phase trajectories for i>( 1 , 2 )(k), k — 2 and 4, of Run A. In this figure, the 
annular phase trajectories (of the coefficients with k x ^ 0) also represent Alfven waves on 
i/ 2 )(0, 1). The mean magnitudes of tA 2 )(k), k = (1, 1) and (-1,1), are significantly larger than 
expected (11.9 vs 1.35), and the mean magnitude of u^(l, 1) is also significantly larger than 
expected (1.74 vs 0.57). Turning to the k 2 = 4 coefficients, he mean magnitudes of v^(k), 
k = (0,2) and (2,0), are also significantly larger than expected (5.39 and 5.78, respectively, 
vs 1.35), while the mean magnitude of 7 /^( 0 , 2) is also significantly larger than expected (1.93 
vs 0.56). 

An example of more ordinary behavior (which absolute equilibrium ensemble theory ex- 
pects to apply all coefficients) is seen in Fig. 7, where the phase trajectory of i)W(k), for 
k = (0, 1) from Run D is given. The mean seems to be around zero and the standard 
deviation of 1.22 appears to be about the predicted value of 1.18. However, if we look at 
i/ 2 )(k), k — 1, anomalous behavior is again observed. In Fig. 8, we see phase trajectories for 
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FIG. 5: Phase trajectories of t>0’ 2 )(k) with k 2 = 1, for Run A. 

t/ 2 )(k), for k = (0, 1) and (1,0), from Run D, with computed mean magnitudes of 8.46 and 
10.30, respectively, differ from the expected mean magnitude of 9.50 (though the rms of the 
computed means is 9.42). 

Our last figure, Fig. 9, draws on data from Run J, whose spectra are shown in Fig. 2. 
These spectra peak at k 2 = K 2 = 36,37, so we look at the phase trajectories of n^(k), 
for k = (39, 46) and (-39,46) as representative examples. There appears to be significant 
structure in Fig. 9, though not as stationary as some of those seen in Figs. 5 and 6 for Run 
A, and Fig. 8 for Run D. In Fig. 9, the coefficients v®(k), k = (36,46) and (-39,46), have 
mean magnitudes 60.4 and 56.1, respectively, while the expectation is 54.0. However, we 
could set our initial conditions so that the u^(k) have energies as small as desired, and in 
this way create essentially static structures in any run, whether their spectra peak at k — 1 
or at k = K. 

We have now seen several examples of dynamical behavior that does not match ensemble 
expectations. We explain this anomalous behavior in the next section, using the concept of 
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— v (2) (-l,l) 
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FIG. 6: Phase trajectories of r/ 1,2 )(k), k 2 = 2 and 4, for Run A. 



FIG. 7: Phase trajectory of uW(k), k = (0, 1), for Run D. 
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FIG. 8: Phase trajectories of the ?/ 2 )(k) with k 2 = 1 for Run D. 



FIG. 9: Phase trajectory of t?( 2 )(k), k = (±39,46), for Run J. 


broken ergodicity. We then discuss the implications for dissipative MHD turbulence. 
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VI. BROKEN ERGODICITY 

The eigenvariables (38) are, again (with C k = \J ft 2 + 7 2 //c 4 and G k = C k — 7/A; 2 ) 

^ (1) ( k ) = ^w(k) + ^fca(k) , (69) 

£ (2) ( k ) = ^2C k G k ~V ^ + Gkk “W • ( 70 ) 

As (45) shows, the expectation values ^|iA 2 )(k)| 2 ^ are inversely proportional to the eigen- 

(n\ (n\ 

values A)y(k), so that when certain eigenvalues A)y(k) are very small, the corresponding 
h (2 )(k) can have very large magnitudes, as Fig. 2 shows, at either lowest or highest k, i.e ., 
at k — 1 or at k — K. 

Whenever (||'fA 2 )(k)| 2 ^ ^(^(k)! 2 ^, equations (69) and (70) can be approximated by 

0 ~ c k G k + ’ ( 71 ) 

6 (2) ( k ) ~ y/2C k G k ^( k ) + G k ka(k) . (72) 

This leads to the approximate relationship G k ip( k) ~ — /?a(k); as \^\/k 2 — > 0, then -0(k) — * 
±a(k) and the corresponding modes no longer interact nonlinearly in the equations of motion 
(1) and (2). 

In general, ^|tA 2 )(k)| 2 ^> ^|i)W(k)| 2 ^ is best satisfied at either k — 1 or k = K and 

at these values of k, ^|f^ 2 ^(k)| 2 ^> is much larger than for any other eigenmode (of the first 
or second kind) at any other values of k. Run A is an example of the largest expected 
eigenmode occurring at k — 1, and Run J is an example of the largest expected eigenmode 
occurring at k — K, as is shown in Fig. 2. In these runs, once the dynamical system comes 

close to equilibrium, and eigenmode energies come close to their expectation values, then 

nonlinear interactions are greatly depressed and the time-evolution of the tA 2 ^(k) for k — 1 
(Run A) or k = K (Run J) slows down to a more or less stationary state, buffeted (slightly 
for Run A and more so for Run J) by what is effectively random noise. This process is seen 
in Fig. 5 for Run A, where fA 2 ^(k) for k = (0, 1) has reached a high degree of stationarity, 
and provides a ‘mean magnetic field’ on which the eigenmode tA 2 ^(k) for k = (1,0) behaves 
like an Alfven wave; further stationarity and Alfven waves are seen in Fig. 6 for some k = 2 
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and 4 eigenmodes in Run A. The process is also seen, to a lesser extent, in Fig. 8 for Run D, 
where the n^(k) for k — 1 are quasi-staionary, and in Fig. 9 for Run J, where the n^(k) 
for k = (±39,46), have reached rougher levels of stationarity. (Here, we have a numerical 
realization of the basic assumption from which Kraichnan started in his prediction of the 
inertial range energy spectrum for MHD turbulence 26 .) 

This process, in which an eigenvalue becomes so small that (71) and (72) are fairly well 
satisfied and one or more eigenmodes become quasi-stationary, so that the expectation of 
zero mean for those eigenmodes is not met in a reasonable time, is the cause of broken 
ergodicity in numerical models of ideal, 2-D, homogeneous MHD turbulence. However, while 
broken ergodicity at k — K is an interesting phenomenon in numerical models of ideal 2-D 
MHD turbulence, if v and rj are nonzero in (1) and (2), then eigenmodes at or near k = K 
will be highly dissipated, obviating the importance of any stationarity. In contrast, though 
dissipation is important at high k, it may have only minimal direct effects at low k for large 
N, in which case we may expect that the ideal results presented here remain pertinent for 
real MHD turbulence. 

Since K may be made larger in 2-D simulations than in 3-D, work is underway to explore 
2-D dissipative MHD turbulence on as large a grid as is practicable, in order to more fully 
investigate broken ergodicity in high Reynolds number MHD turbulence. Nevertheless, it is 
also important to continue 3-D studies, as there are quantitative and qualitative differences 
between 2-D and 3-D MHD turbulence. Ideal 3-D MHD turbulence exhibits broken ergodicity 
at k — 1, but not at k — 2 or 4, though this may be due to the fact that long-term 
simulations 10 have used 32 3 grids, with K = 15.08, while in the 2-D simulations presented 
here, broken ergodicity appears at k — 1, 2, 4, and higher, when (p 0 > §, but here 128 2 
grids were used, with K = 60.34. Thus, still larger grids, which may more fully resolve the 
details of intrinsic coherent structure, are necessary for both 2-D and 3-D MHD turbulence 
simulations, remembering that very long simulation times are also required. The challenge 
is to find the most efficacious combination of constantly evolving computer software and 
hardware available for the task. 

A major qualitative difference between ideal 2-D and 3-D MHD turbulence is that ideal 
spectra will peak at k — K for (p Q < ± 2-D runs, while this cannot happen in ideal 3-D MHD 
simulations (expectation values of magnetic energy are never less than those for kinetic 
energy in ideal 3-D MHD turbulence 28 , but more pertinently, the smallest eigenvalues occur 
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at only at k — 1 for any expected value of magnetic energy 10 ). However, there are enough 
similarities between 2-D and 3-D MHD turbulence that both merit continued study, both 
for their inherent mathematical interest and for their important applications. 


VII. DISCUSSION 

Broken ergodicity in 2-D MHD turbulence depends on the approximations (71) and (72) 
being good ones that produce the relationship Gkip( k) ~ —/3a( k) for either k = 1 or k = 
K. This relationship between ^(k) and a(k) involves the cofactors Gk and f3 and thus 
depends explicitly on the inverse temperatures f3 and 7 . In contrast, analogous results for 
3-D ideal MHD turbulence 10 lead to relationships (but only at k — 1) such as u(k) ~ a>(k), 
j(k) « b(k) and u(k) ~ b(k), where |k| = 1 , that do not contain cofactors involving inverse 
temperatures. However, these relationships are similar in that they indicate that both 2-D 
and 3-D magnetofluids relax to ‘force-free’ states. When dissipation is present, higher k 
modes decay quicker, so most of the energy eventually winds up in the k — 1 modes, which, 
if these ideal results are pertinent, become force-free states of the maximum size permitted 
by boundary conditions. 

Let us briefly place the matter into a historical context. In 1974, Taylor put forth the 
idea that a cylindrically confined plasma somehow relaxed to a state in which energy E was 
minimized while magnetic hclicity Hm remained constant 2 '. This led to a force- free config- 
uration described by lowest-order Bessel functions. The qualitative process through which 
these Taylor states occur has been termed ‘selective decay ’ 29 or ‘dynamic alignment’ 30 , both 
of which were seen later to be synonymous with ‘plasma relaxation’ 31 . Related concepts 
are ‘self-organization’ in a plasma 32 and ‘depression of nonlinearity’ 33,34 . ‘Broken ergodic- 
ity’, in turn, provides a quantitative, mathematical explanation for the process of large-scale 
structure formation, dynamo action and self-organization in model systems describing ideal, 
homogeneous MHD turbulence. The addition of dissipation to these model systems has been 
seen through numerical simulation 16,28 to lead to selective decay and dynamical alignment, 
while the quasi-stationary large-scale structures due to broken ergodicity were still in evi- 
dence. Although these simulations were on only moderately sized grids, we fully expect to 
see even more compelling evidence as future long-time simulations that use larger and larger 
sized grids become possible. 
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VIII. CONCLUSION 


In summary, we have explained previous observations of coherent structure and apparent 
non-ergodicity in numerical simulations of 2-D MHD turbulence as being due to broken 
ergodicity. As stated earlier, the next step is to perform dissipative simulations on as large 
a grid as practicable, to see how ideal results translate over to more realistic model systems. 
Keeping these simulations in a 2-D framework allows us to keep the maximum wavenumber 
K, as well as the dissipation wavenumber ku < K, as large as possible. The results presented 
here indicate that there is a rich and complex nonlinear dynamics inherent in the model 
systems that represent 2-D MHD turbulence, and further study should reveal even more 
detail. 

Advances in 3-D simulations are, perhaps, of more practical importance as these sim- 
ulations are more directly related to understanding turbulent plasmas in both laboratory 
and astrophysical settings. However, since maximum wavenumbers K are always going to 
be greater in 2-D simulations, the higher Reynolds numbers that are possible will continue 
to make 2-D simulations attractive for the study of MHD turbulence, notwithstanding the 
inherent difference between the 2-D and 3-D cases. The increasing capability of computer 
systems with regard to size and speed will allow us to gain even more insight into both 2-D 
and 3-D MHD turbulence, and we look forward to learning of new results in both cases. 
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